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We develop new variational principles to study the structure and the stability of equilibrium states 
of axisymmetric flows. We show that the axisymmetric Euler equations for inviscid flows admit 
$H . an infinite number of steady state solutions. We find their general form and provide analytical 

solutions in some special cases. The system can be trapped in one of these steady states as a 
result of an inviscid violent relaxation. We show that the stable steady states maximize a (non- 
universal) //-function while conserving energy, helicity, circulation and angular momentum (robust 
constraints). This can be viewed as a form of generalized selective decay principle. We derive 
relaxation equations which can be used as numerical algorithm to construct nonlinearly dynamically 
stable stationary solutions of axisymmetric flows. We also develop a thermodynamical approach 
to predict the equilibrium state at some fixed coarse-grained scale. We show that the resulting 
distribution can be divided in two parts: one universal coming from the conservation of robust 
I ' invariants and one non-universal determined by the initial conditions through the fragile invariants 

(for freely evolving systems) or by a prior distribution encoding non-ideal effects such as viscosity, 
small-scale forcing and dissipation (for forced systems). Finally, we derive a parameterization of 
^ ' inviscid mixing to describe the dynamics of the system at the coarse-grained scale. A conceptual 

, interest of this axisymmetric model is to be intermediate between 2D and 3D turbulence. 
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The ubiquity of rotating systems in astrophysics and geophysics makes axisymmetric flows a classical paradigm. 
In the laboratory, two axisymmetric devices, the Taylor-Couette flow and the von Karman flow, have become a 
standard to investigate issues such as super-critical and sub-critical stability [| , fluctuation of global quantities 0, Q 
or turbulent transport However, many basic issues regarding stability and turbulence in axisymmetric flows 

still remain unsolved. For example, one still fails to understand the onset of turbulence or the equilibrium state 

• i-h | in Taylor-Couette with outer rotating cylinder Q, or the recent bifurcation of the turbulent state observed in von 

Karman flow Q . 

In the past, dynamical stability and equilibrium properties of flows have often been studied using variational [9j or 
fS maximization [lOj principles. Examples of application to axisymmetric flows include necessary criteria for instability or 
turbulent velocity profiles in Taylor-Couette flow. Maximization or minimization principles have also been used to give 

• *-h [ sufficient criteria of nonlinear dynamical stability One interest of these methods is their robustness, in the sense 

that they mostly depend on characteristic global quantities of the system (such as the energy) but not necessarily on 
small-scale dissipation or boundary conditions. More recently, optimization methods have been developed within the 
framework of statistical mechanics for two-dimensional (2D) perfect fluids. In that case, variational principles based 
on entropy maximization determine conditions of thermodynamical stability. Onsager [12j first used a Hamiltonian 
model of point vortices and identified turbulence as a state of negative t emp erature leading to the coalescence of 
vortices of same sign [l^. Further improvements were provided by Kuzmin |l4j . Miller 15] and Robert and Sommeria 
|16| who independently introduced a discretization of the vorticity in a certain number of levels to account for the 
continuous nature of vorticity. Using the maximum entropy formalism of statistical mechanics 17], it is then possible 
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to obtain the shape of the metaequilibrium solution of Euler's equation as well as the distribution of the fine-grained 
fluctuations around it. A variety of solutions are found and the bifurcation diagram displays a rich structure as 
illustrated by Chavanis & Sommeria [I^ m a particular limit of the statistical theory. Two-dimensional turbulence is 
however very peculiar since it misses vortex stretching, one essential ingredient in 3D turbulence. The adaptation of 
these methods to more realistic situations is therefore not obvious. 

In the case where the system admits a scalar invariant (Dt<J = dt<J + u • Vct = 0), one can show that a Liouville 
theorem holds (incompressibility of the motion in phase space). Indeed, the proof given by Kraichnan & Montgomery 
|19| in the case of 2D turbulence can be extended to any dimensional turbulence with a conserved quantity. This is 
in fact the case for axisymmetric flows where the symmetry imposes angular momentum conservation. Due to violent 
relaxation, the system is expected to reach a metaequilibrium state which is a steady solution of the axisymmetric 
Euler equations. The purpose of the present paper is to derive explicit results regarding nonlinear dynamical stability 
and thermodynamical equilibrium properties of axisymmetric flows using optimization methods. 

In the first part of the paper (Sec. we consider the nonlinear dynamical stability of stationary solutions of 
the axisymmetric Euler equations (i.e. without viscosity). These equations are written in Sec. Ill Al and the general 
form of stationary solutions is obtained in Sec. Ill Bl In Sec. Ill CI we list the conservation laws of the axisymmetric 
Euler equations. We find non trivial invariants in addition to the usual ones. In Sec. Ill Dl we show that the 
equilibrium solutions can be obtained by extremizing a functional built with all the invariants. The fact that this 
optimization procedure returns the general stationary solution means that we have found all the invariants. In 
Sec. Ill El we distinguish between fragile (Casimirs) and robust (energy, helicity, circulation, angular momentum...) 
invariants. We argue that, in the presence of viscosity or coarse-graining, the metaequilibrium state maximizes a 
certain (non-universal) H-function while conserving the robust constraints. This is similar to the case of pure 2D 
hydrodynamics pfl . l2lll22ll23j | except for the replacement of vorticity by angular momentum. In Sec. Ill Fl we propose 
a numerical algorithm based on the maximization of the production of an H-function while conserving the robust 
constraints. This can be used to compute numerically arbitrary nonlinearly dynamically stable stationary solutions 
of the axisymmetric Euler equations. This is similar to the relaxation equations proposed by Chavanis [23, ISl 
in pure 2D hydrodynamics. In Sec. Ill Gl we provide simple analytical steady solutions of axisymmetric equilibrium 
flows. In the second part of the paper (Sec. IIII|I . we develop the statistical mechanics of such flows to predict the 
metaequilibrium state. The statistical equilibrium state is obtained by maximizing a mixing entropy (Sec. IIII A|) 
while taking into account all the constraints of the dynamics. This yields a Gibbs state (Sec. 1111 Bfl which gives the 
equilibrium coarse-grained angular momentum as well as the fluctuations around it. We check that the coarse-grained 
field is a stationary solution of the axisymmetric Euler equations. However, since the Casimirs are not conserved on 
the coarse-grained scale, the distribution of fluctuations is non universal and depends on the initial conditions (or 
fine-grained constraints). This is also the case for the coarse-grained field. In Sec. IIII Dl we use a maximization 
of the entropy production to derive relaxation equations towards the statistical equilibrium state. This is similar to 



the approach proposed by Robert & Sommeria | 24| in pure 2D hydrodynamics. Finally, in Sec. IIII El we introduce 
the notion of prior distribution of fluctuations for systems that are forced at small scales. We show that the coarse- 
grained field maximizes at statistical equilibrium a generalized entropy fixed by the prior distribution. Thus, the 
relaxation equations based on the maximization of the production of generalized entropy (Sec. Ill F|> can also provide 
a parameterization of axisymmetric turbulence in the presence of a small-scale forcing |20t l2ll |22|. 

II. DYNAMICAL STABILITY OF AXISYMMETRIC FLOWS 

A. The axisymmetric Euler equations 

The Euler equations describing the dynamics of an inviscid incompressible axisymmetric flow can be written in 
cylindrical coordinates (r, 9, z) as: 

-d r (ru)+d z w = 0, (1) 
r 

v 2 1 

dtu + ud r u + wd z u = — d r p , 

r p 

dtv + ud r v + wd z v H = , 

r 

dtw + ud r w + wd z w = d z p , 

P 

where (u,v,w) denote the components of the velocity in a cylindrical referential. Note that the third equation 
expresses the conservation of the angular momentum a = rv. The two other equations for u and w involve a pressure 
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field determined through incompressibility. However, it can be eliminated by using the stream-function vorticity 
formulation |25j. The two new scalar variables are the azimuthal component of the vorticity ujg = d z u — d r w and the 
stream function ip defined by: 

u = — d z ip, and w — —d r i\) . 
r r 

The existence of a stream function results from the incompressibility and the axisymmetry of the flow. In this 
formulation, the system Q can be rewritten: 

d t v - ^-d z ipd r v + ^d r ipd z v - -^dztfjv — , (2) 

1 1 LOg ( v 2 

d t ui d z ipd r u! 6 + -d r ipd z uj g + d z i>— = d z \ — 

By definition, the azimuthal component of the vorticity is related to the stream function by: 

d r (~d r A + -%1> = -UJg . (3) 

\r J r 

We now introduce two new fields, the angular momentum a — rv and £ which is related to the azimuthal component 
of the vorticity by £ = ojg/r. Changing variables from (r, z) to (y, z) where y — r 2 /2, we can finally recast the equations 
© and J3J as 



d t <j + {yj,<j} = 0, (4) 

Ay 2 



dtt + Mt} = d z [— 2 



where {ip, 0} = d y i\)d z <i> — d z ipd y (j) is the Jacobian and A* is a pseudo-Laplacian. We also note that u z = d y tp = w 
and u y — —d z ip — ru. This formulation of the axisymmetric Navier-Stokes equation has to be supplemented by 
appropriate boundary conditions. For reasons which will become clear later, we delay this topic until the discussion 
of the conservation laws (section III CI and Appendix ^ . In the following, we study axisymmetric equilibrium flows 
by using the system of equations |@J instead of (JIJ . We look for stationary solutions and investigate their stability 
by using variational methods. Notice that only two scalar variables are sufficient to prescribe such flows: we use a, 
related to the azimuthal component of the velocity field and £, related to the azimuthal component of the vorticity. 

B. Stationary solutions 

We now derive the general form of stationary solutions of the axisymmetric Euler equations (@J) . Noting that 

o , „ fa 2 

I 

the stationary equations can be written 



{M = and {V, £} + <j} = 0. (6) 

The first equation is satisfied if 

V; = R(a), (7) 
where R is an arbitrary function. Using the general identity 

{R(a), a - R'(a){a, £} = {a, £&(*)}, (8) 

the second equation becomes 

{a^R'(a)} + {^- ,a} = 0, (9) 
22/ 
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or, equivalently 



{*,&(*) -£} = 0. (10) 



Therefore, the general stationary solution of Eqs. 10} is of the form 

ip = R(a) and £R'(a) - = G(a), (11) 

where R and G are arbitrary functions. If R is monotonic, we can set / = R^ 1 and we get a = f(ip) and 

*--2f R>[R-W)] =9W - (1) 

Using the identity 
we finally obtain 

_ A ^ = ^ = M/' W+5W) (14) 

where / and g are arbitrary functions. We can obtain these equations directly if we note that Eq. (JSJ-a) is satisfied if 
a = f(ip)- Then, using the general identity 

we can rewrite Eq. JSJ-b) in the form 



22/ J I 2y 



tf-Sffl}^., (16) 

which leads to Eq. H14l b1 . Equation Ijl^ b) is the fundamental differential equation of the problem which must be 
supplemented by appropriate boundary conditions. Some particular solutions of this equation will be given in Sec. 
Ill Gl We will first show that the stationary solutions can be found by a variational principle depending only on the 
conservation laws of the system. 



C. Conservation laws 



Axisymmetric inviscid flows satisfy a number of conservation laws. We here give the expression of these conserved 
quantities and postpone corresponding proofs in Appendix^ To derive the two first conservation laws, we must 
assume that the function if) vanishes on the boundary of the domain which amounts to considering that the normal 
component of the velocity is zero at the boundary. This condition is not sufficient for deriving the third conservation 
law and one must also suppose that either a or £ vanishes at the boundary. 



The first conserved quantity is the total energy 



E = i J(u 2 +w 2 )rdrdz + ^ J v 2 rdrdz (17) 
= \ J ^ d V dz + \ J ~ d v dz - 



Here, we have normalized the energy by 2-7T and used integration by parts to obtain the second expression. 



5 



• Because of l|H-al any function of the angular momentum is also an invariant noted as 

If = J f (a) dydz. (18) 

These functionals are called the Casimirs. The conservation of all the Casimirs is equivalent to the conservation 
of all the moments of cr, denoted I n — J a 11 dydz. The first moment is the total angular momentum I = J a dydz. 
If a = or, more generally, a = a(y), then £ is conserved via (Q}b). In that case, £ is called a potential vorticity 
(or a pseudo-vorticity) and there is an additional class of Casimir invariants: Ih = J dydz. We ignore this 
difficulty linked to a sort of "degeneracy" for the time being. Note that the situation where only the pseudo- 
vorticity is conserved (i.e the case a = 0) has been treated in [2|j. In that case, the generalization essentially 
amounts to replacing the Laplacian A in pure 2D flows by the pseudo-Laplacian A*. The situation that we 
consider here is complicated by the existence of additional invariants such as helicity discussed below. This 
makes our situation intermediate between pure 2D turbulence and 3D turbulence, an interesting feature of our 
model. 

• Like in any 3D flows, the total helicity, H = Jv ■ ujrdrdz = J dydz, is also an invariant. However, more 
generally, we show in Appendix [S] the conservation of a generalized helicity 

H F = J £F(a) dydz , (19) 

where F is an arbitrary function. In particular, the total vorticity Hq = T = J £ dydz is conserved. 



D. Nonlinear dynamical stability 

From the integral constraints discussed previously, we can build a functional F = E + If + Hp. This functional 
is an invariant of the inviscid dynamics. This is an extension of the Energy-Casimir functional considered in |llj| . It 
is also similar to a free energy in thermodynamics. We now show that a critical point of F determines a stationary 
solution of the axisymmetric Euler equations. Furthermore, following [Tl| . a minimum or a maximum of F provides 
a condition of formal nonlinear dynamical stability. This means that a perturbation will remain close (in some norm) 
to this minimum or maximum. Writing 

SF = 5{E + I f + H F ) = 0, (20) 

and taking variations on a and £, we obtain 

i/j + F(a) = and + f(a) + £F'(a) = . (21) 
22/ 

Setting R = —F and G — f , we recover the equations (|llf) characterizing a steady solution of the axisymmetric Euler 
equations. Since we obtain the general form of steady states it means that we have found all the conservation laws of 
the axisymmetric Euler equations. 

In order to gain some physical insight in the problem, we consider from now on a simpler model where only the 
usual helicity H and the total vorticity T are conserved instead of all the generalized helicities. This is similar to our 
choice of restricting ourselves to the Chandrasekhar model in axisymmetric MHD (2?J . We define 

S[tr] = - J C{a) dydz , (22) 

where C is an arbitrary convex function, i.e. C" > 0. Such functionals are exactly conserved by the axisymmetric 
equations (they are particular Casimirs). Therefore, as in 2D hydrodynamics |28| . the maximization of S at fixed 
energy E, helicity H, circulation T and angular momentum / determines a nonlinearly dynamically stable stationary 
solution of the axisymmetric Euler equations. This refined stability criterion is stronger than the maximization of 
J = S — f3E — fiH — vT — al which just provides a sufficient condition of formal nonlinear dynamical stability |ll| . 
The difference between these two criteria is similar to a notion of ensemble inequivalence in thermodynamics (where 
S plays the role of an entropy and J the role of a free energy) |2(3, . We shall not prove the nonlinear dynamical 
stability result in this paper and refer to [2^ for a precise discussion in 2D hydrodynamics. In Sec. Ill El we show, 



6 



however, that this maximization principle is consistent with the phenomenology of axisymmetric turbulence provided 
that a is interpreted as the coarse-grained angular momentum. 
To first order, the variational problem takes the form: 

SS - f35E - [iSH - 7 <5r - aSI = 0, (23) 

where a, (3, fi and 7 are appropriate Lagrange multipliers. This variational problem determining nonlinearly dynam- 
ically stable stationary solutions of the Euler equations is similar to a variational problem in thermodynamics where 
S plays the role of an entropy and (5 the role of an inverse temperature Using the expression of S, E, H, T 

and /, we find that the solutions of ((23(1 valid for any da and 6£ satisfy 

M = -iia - 7, (24) 

which is a particular case of Eq. Thus, the variational principle selects stationary solutions of the axisymmetric 

Euler equations. We note that when only the ordinary helicity is considered (instead of the general helicity), we obtain 
a linear relationship between a and ip. This is similar to the linear relationship between velocity V and magnetic field 
B in MHD |27j . Note that we have just considered the first order variations here. To check if solutions 1)24(1 are real 
maxima of S, one has to look for second-order variations as discussed in Appendix iBl 



E. //-functions and generalized selective decay principle 



We now introduce the notion of fine-grained and coarse-grained fields. The first one refers to the original field 
defined on all points of space and time and the second one to a smooth version of it, where a local average of the 
field has been performed. The coarse-grained field is also defined on every point of space but contains less small-scale 
structure than the original field. Since the functionals l|22f) calculated with the fine-grained field a are particular 
Casimirs, they are rigorously conserved by the fine-grained dynamics. In contrast, as Tremaine et al. [30| have shown 
for the Vlasov equation in stellar dynamics, the functionals of the form 1(22(1 calculated with the coarse-grained field 
increase as mixing proceeds. This is similar to the Boltzmann //-theorem in kinetic theory except that the Vlasov 
equation does not single out a particular functional. This property is true also in the present context since Eq. 
(0}a) plays the same role as the Vlasov equation. Therefore, S[a] = — J C(a)dr increases along the dynamics in the 
sense that S[a(r,t)] > S\a(r,0)] for all C and all t > where it is assumed that, initially, the flow is not mixed: 
cx(r,0) = a(r 0) (note that nothing is implied concerning the relative values of S(t) and S(t') for t > 0, t' > 0). 
Following [3(j, these functionals will be called //-functions (or generalized //-functions) . They also increase (in that 
case monotonically) in the presence of viscosity since the equations of motion now become: 



d t a + {i/j,a} = vr \ - ) = /' 



2 da 
Aa-- — 

r or 



(25) 



a 



r dr 



and by integration by part, one can show that S — v J C"(a)(Va) 2 dv > 0. By contrast, the integrals E, T, H, and 
/ are exactly or approximately conserved on the coarse-grained scale (i.e. when they are calculated with the coarse- 
grained field) or in the presence of a small viscosity. For example, in the presence of viscosity, the kinetic energy evolves 
such that E = —v J u> 2 dr. It is easy to show that for axisymmetric fields, the total vorticity uj — r£eg + V x (a/reg) 
vanishes in the long time limit. The demonstration is similar to Cowling's [3lj theorem of dynamo theory which 
states that an axisymmetric magnetic field cannot grow in an axisymmetric velocity field: the first equation in 1251) 
shows that a — * for large time and, consequently, the source term in the second equation d z (<t 2 /4j/ 2 ) vanishes 
in the long time limit, which implies that £ — > 0. Thus, for axisymmetric flows, both components of the vorticity 
vanish in the long time limit and the energy is approximately conserved. In a similar way, it can be shown that T, H, 
and / are approximately conserved and must therefore be strictly taken into account in the constraints. Therefore, 
the functionals S can be viewed as fragile invariants while the constraints E, T, H, and I are robust invariants. 
This generalizes the notion of selective decay in pure 2D turbulence where the enstrophy decays while the energy is 
approximately conserved. In fact, minus the enstrophy is a particular H- function |2Ct l2ll |23|. The same discussion 
applies in the present context. On the basis of this phenomenological principle, we expect that, due to chaotic mixing 
and violent relaxation, the system will reach a metaequilibrium state which maximizes a certain //-function (non- 
universal) at fixed E, T, H, and /. This phenomenological argument returns the variational principle l(23|) . Since 
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this metaequilibrium state results from turbulent mixing, it is expected to be particularly robust and should possess 
properties of nonlinear dynamical stability. Therefore, the stability arguments given previously are consistent with the 
phenomenology of axisymmetric turbulence, provided that a is interpreted as the coarse-grained angular momentum 
W. This is remarkable because the two arguments are relatively independent: there is no direct notion of decay (of 
—S) in the first argument while this lies at the heart of the second. In fact, the phenomenology of violent relaxation 
explains how an inviscid system can reach a nonlinearly dynamically stable stationary state on a coarse-grained scale 
which is a maximum of a certain TJ-function at fixed robust constraints (while S[a] is rigorously conserved on the 
fine-grained scale). The point is that during mixing Da/ Dt ^ and the _ff-functions S\a] increase. Once it has mixed 
Da/ Dt = so that S\a] = 0. If a(r, t) has been brought to a maximum efo(r) of a certain H- function (as a result of 
mixing) and since S[a] is conserved (after mixing), then ef is a nonlinearly dynamically stable stationary solution of 
the axisymmetric Euler equation according to the stability criterion of Sec. Ill Dl 



F. A numerical algorithm for the dynamical stability problem 



We shall construct a set of relaxation equations that increase S[a] while conserving all the robust constraints E, 
r, H, and /. These relaxation equations, which solve the optimization problem of Sec. Ill Dl can serve as powerful 
numerical algorithm [44j to compute arbitrary stationary solutions of the axisymmetric Euler equations. In addition, 
they guarantee that these solutions are nonlinearly dynamically stable with respect to the inviscid dynamics. Such 
relaxation equations therefore have a clear practical interest. They extend those obtained by Chavanis pel I2H l2^ | in 
2D hydrodynamics. 

We write the dynamical equations as 

! = -v-J„, § = -v-Je, (26) 

where 3 a and are the currents to be determined. We have not added advective terms since we here use these 
equations just as numerical algorithms, not as a description of the dynamics (see, however, Sec. IIII Efl . By construction, 
these equations satisfy the conservation of the total vorticity and total angular momentum. On the other hand, the 
conservation of energy and helicity impose the constraints 

1 



E 



= JfVipdr 



H = = 



V£ dr 



• Vcr dr. 



Finally, the time variations of S[a] are given by 

/ C"{a) J a ■ Vcr dr. 



We derive the optimal currents which maximize S with E = H = and the additional constraints 

^<C CT (r,i) and ^<Q(r,t), 



putting an upper bound on the currents. Writing the variational principle in the form 

SS - f3(t)SE - fi(t)8H - 
we obtain the optimal currents 



\ H y \dr 



C"(a) 



(27) 



(28) 



(29) 



(30) 



(31) 



(32) 



where (3{t) and /i(t) are Lagrange multipliers which evolve in time so as to conserve energy and helicity. They are 
determined by substituting l|32|l in the constraints (|27Jl and (|28|) . Plugging the optimal currents in H26JI we get 



dt 



Vcr 



(3 



C"(a) \2yJ C"(a) 
D'(/3V^ + A*Vcr) 



(33) 
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It is straightforward to check that the system (|33|l increases the functional 1|22|) . i.e. S = J (J 2 /Dp + J^/D'jdr > 0, 
and that the stationary state is given by Eqs. i|24|) . Using the same principle, we can also write relaxation equations 
which minimize T . The optimal currents are 



W + ^) V fe) + ^) V( ^' M) 



(34) 



J, = -D 
J c = -D'(Vip + VF(a)). 
and they return as an equilibrium state, the stationary solutions (|21H . 

G. Analytical solutions in simple cases 

The steady state equations <|14|) admit analytical solutions for simple shapes of the arbitrary functions / and g. We 
will here derive some of these solutions and show that they are critical point of simple functionals. 

1. g = 

Let us first consider the case where g = 0. In that case, the steady solution obeys 

* = fU>) and - A*V = £ = ^P-f'm • (35) 
This equation admits simple solutions independent of y. Indeed, the second equation becomes 

d^ = "2^ (/) ' (36) 

which is equivalent to the motion of a particle in a potential ^f 2 (ip) where ip plays the role of position and z the 
role of time. Multiplying both sides of Eq. I|36|l by dip/dz, and then integrating twice, the solution can be put under 
parametric representation as 

J \/A'i - P{4.) 

where K is an integration constant and we have returned to original variables. For example, for linear /, one obtains 
a oc ip oc cos(Kz). 

2. Constant g and linear f 

Consider now the case where g is a constant g = C and / is a linear function of ip, f = A + Bip. The equations 
become 

a = A + BtP , (38) 

AB + Bhp 
-A,V = £ = h C . 

Note that these equations arise as critical points of the functional J-q = E + [IqH + vqT + uqI, i.e. they determine 
a state of minimum energy at fixed H , T and /. Equation (|38l -b) is an inhomogeneous linear equation for ip. The 
general solution is the sum of a special solution of the inhomogeneous equation superposed to the general solution of 
the homogeneous equation. A special solution is easily found as 
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This solution corresponds either to a laminar Taylor-Couette profile, or to a profile maximizing turbulent transport 
in Taylor-Couette flow 01 ■ Notice that the present theory is unable to capture the 1/r dependence of the Taylor- 
Couette flow because the solutions have to be regular at the origin. To reproduce such a behavior, one has to consider 
a domain with two boundaries (corresponding to the inner and outer cylinders), one with ip = and the other with 
ip 0. Such a procedure would introduce boundary terms in the conserved quantities. 

A general solution of the homogeneous equation can be found by the method of separation of variables by writing 
ip = G{y)F(z). It is then easy to show that F(z) oc cos(kz + <p) where k and <fi are two constants. Then, one finds 
that G obeys the following equation: 



d 2 G 
dy 2 



B 2 — n A 

2y 



G = 



dr 2 



dr 



+ [(B 2 - n 2 )r 2 



1 -=0 

r 



(40) 



whose solution can be expressed in term of Bessel function of first order. The general solution for ip is thus 



jirJi (y/ [B 2 — n 2 )rj cos(kz + . 



(41) 



where n, k and are integration constants. This solution is a critical point of the functional Tb = E + /j,qH, i.e. it 
minimizes the energy at fixed Helicity. In vectorial form, this leads to curl(u) = Au such that vorticity and velocity 
are everywhere proportional. This is the so-called Beltrami solution which has proven to be important in the study 
of the dynamo effect [12, ie. the generation of a magnetic field by a conducting fluid. The most popular flow of this 
type is the Roberts flow [23|. In the limit k — > B, the homogeneous solution tends to r 2 and this solution becomes 
equivalent to one of the nonlinear self-similar solution of the von Karman flow found by Zandbergen and Dijkstra 
[34|. On figure ^ we show a contour plot of this solution. 



V (r.z) 





FIG. 1: (Color online) Contour plot of the stream function associated to the "Beltrami solution"for the parameters A — 3, 
B — 2, C — 2, n — 1/2, k, — ir, and 4> — 0. The left hand side corresponds to the homogeneous solution 14111 and the right hand 
side to the sum of the homogeneous solution and the particular solution 13911 . 



3. f and g linear 

The case where both / and g are linear g — C + Dip and / = A + Bip is similar to the previous one. The equations 
are now 



a = A + Bip, 
' B 2 



[2y 



D 



(42) 
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Note that these equations are obtained by minimizing the second moment of angular momentum I2 = J a 2 dydz at fixed 
E, H, r and /. They represent, therefore, the counterpart of the minimum enstrophy principle in 2D hydrodynamics, 
leading to a linear relationship between vorticity and stream function. Equation H42l b) is an inhomogeneous linear 
equation for ip. The general solution is the sum of a special solution of the inhomogeneous equation superposed to 
the general solution of the homogeneous equation. Solutions of the homogeneous equation can be found by assuming 
separation of variable as previously t/j = F(z)G(y). The solution for F is F{z) = cos(kz + <fi), where k and <j> are two 
integration constants. Then, the equation for G involves a supplementary term compared to the one in the previous 
section: 



d 2 G 
dy 2 



B 2 



2y 



D] G = 0. 



equations 

-2' 



can be expressed in terms of 



Whittaker 

A 



The two solutions of this 

W\ t ±i/2(2^—Dy), A = (B 2 — n 2 )/8D. These function behave at infinity like y x exp( 
variable, one can therefore express the general solution of equation (|42|l as 

■0 O = fiWx,±i/ 2 {V^Dr 2 ) cos(kz + 4>), 



(43) 

function (see |35|. p. 1059) 
-y/2). Turning back to original 

(44) 



where C is an integration constant. Note that the negativeness of the coefficient D is imposed by an asymptotical 
analysis of equation 1)43(1 : when y — > 00, we see that a positive coefficient D would introduce an oscillatory, unphysical, 
behavior for G. Figure shows a typical realization of this solution. 




FIG. 2: (Color online) Contour plot of the stream function associated to the "Whittaker solution" of equation 1431 for the 
parameters A — 3, B = 2, C = 2, D = —3, /U = 1/2, k — n and (f> = 0. 



III. STATISTICAL MECHANICS OF AXISYMMETRIC FLOWS 



A. The mixing entropy 



Starting with a given initial condition, equations Q are expected to develop a complicated mixing process, with 
formation of finer and finer structure, leading to more and more degrees of freedom. A precise prediction of the state 
of the system would a priori require to keep track of all these degrees of freedom. Suppose however that we are only 
interested in the knowledge of the system at some coarse-grained scale. Since mixing is continuously occurring at 
smaller and smaller scales, we can expect the formation of a metaequilibrium state on the coarse-grained scale. Our 
goal is to derive its shape through thermodynamical arguments, based on entropy maximization in the spirit of |17| . 
We focus here on basics. More discussion about this procedure can be found in |l6(. We introduce coarse-grained and 
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fine-grained quantities. According to Eq. Q}c), ip — A^ 1 ^ is expressed as an integral over £. It is therefore a smooth 
function whose fluctuations can be neglected ip — ip. We shall determine the distribution of fluctuations of angular 
momentum a by an approach similar to that developed in 2D turbulence. We then introduce p(r, rf) , the density 
probability of finding the value a = 77 at position r. Then, the coarse-grained angular momentum is a = J pr\dr\ and 
the local normalization J pdr] = 1. We introduce the mixing entropy 

S[p] = — J plnpdydzdrj. (45) 

This functional S[p] is equal to the logarithm of the disorder where the disorder is the number of microstates cor- 
responding to the macrostate p{v,rf) which can be obtained by a combinatorial analysis (see ,2l| for details). The 
fact that the entropy depends only on the distribution of angular momentum guarantees that it will be conserved 
by the fine-grained dynamics (however, it increases on a coarse-grained scale as we will show subsequently). As a 
drawback, we will not be able to characterize the fluctuations of pseudo-vorticity £ by this method; we will simply get 
its coarse-grained value. This difficulty may reflect the fact that we are in a situation intermediate between 2D and 
3D turbulence. In 2D turbulence, the distribution of vorticity is enough to construct equilibrium solutions whereas, 
in 3D turbulence, it is well-known that the average of fluctuating quantities (such as the Reynolds stress) are very 
important. In Appendix [D] we consider another approach which puts the fluctuations of a and £ on an equal footing. 
This approach seems to indicate that the fluctuations of £ have a peculiar behavior that may give rise to a sort of 
"phase transition" . Here, for simplicity and clarity, we first concentrate on the simplified situation (similar to 2D 
turbulence) where the fluctuations of £ are mild; however, we keep in mind that there may be another regime (closer 
to 3D turbulence) characterizing axisymmetric flows. 

The coarse-grained values of the constraints can be written as 

W.,,. I fa 2 



E = - J tydydz + - J — dydz , (46) 



H = Jfadydz, T = j£dydz, (47) 
T n = [ prf dydzd V . (48) 



For the coarse-grained energy, we made the non trivial hypothesis that the fluctuations of the energy could be neglected 
and used the coarse-grained field a to calculate the mean energy. This is justified by the remark made in section III El 
that in the presence of a small viscosity (or a coarse-graining), the energy is approximately conserved. 

As indicated previously, E, T, H and / are robust constraints. Thus they can be determined at any time since their 
exact value is close to their coarse-grained value. They have thus been expressed directly in terms of the coarse-grained 
field. By contrast the higher moments I n >i are fragile constraints because they are altered by coarse-graining. They 
can be determined only from the initial conditions which are supposed un-mixed (or from the fine-grained field) since 
they are affected by coarse-graining in the sense that I n> i\a] = ja^dydz ^ / n >i[cr] = J <j n dydz = J prfdydzdr). 
Part of the Casimirs goes into the coarse-grained field and part goes into fine-grained fluctuations. In a sense, I n >i 
are "hidden constraints" because they cannot be determined from the coarse-grained field. Therefore, the statistical 
theory assumes that we know the initial conditions in detail and that these initial conditions represent the fine- 
grained field. This poses a conceptual problem because in real situations we are not sure whether the initial condition 
is already mixed or not. Furthermore, if the initial condition already results from a mixing process (like vortices 
formed in a succession of mergings in 2D decaying turbulence), it is more proper to ignore its fine structure and 
take its coarse-grained density as a new "fine-grained" initial condition to predict the next merging (see [H, p 284). 
In fact, due to viscosity, the fine structure of the field is progressively erased. These difficulties are intrinsic to the 
statistical theory of continuous fields. 

B. The Gibbs state 

The most probable distribution at metaequilibrium is obtained by maximizing the mixing entropy S[p] at fixed E, 
H, r, /, I n and normalization. We introduce Lagrange multipliers and write the variational principle in the form 

SS - (35E - pSH - -fST -aS7-^T a n S ( f prf 1 dydzdrA (49) 

n>l ^ ' 

- j Ciy,z)s(^J pd-q^j dydz = . 
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The last term in this equation corresponds to the normalization of the probability density in each point of space and 
thus needs the introduction of one Lagrange multiplier C,{y,z) for each point (y,z). We shall treat the variations on 
p and £ independently. The variations on £ imply 



pip = -fia - 7. (50) 



The variations on p yield the Gibbs state 



where 



P = VT—^)e- {f ^ + ^ +a)ri , (51) 
Z(y,z) 



x (rj) = e -^n>x^v n _ (52) 



To prepare the approach of Sec. IIII El we have distinguished between the Lagrange multipliers a n> \ which account 
for the conservation of the fragile constraints /„>i (they have been regrouped in the function x(v)) horn the Lagrange 
multipliers /?, p and a which are related to the robust constraints. Therefore, the distribution j|ST|) is the product of a 
universal Boltzmann factor and of a non-universal function xiv) depending on the initial conditions. Note that instead 
of conserving the fine-grained moments I n = J a n dydz — J prf 1 drjdydz of angular momentum, we can equivalently 
conserve the total area 7(77) = / pdydz of each level. The "partition function" is determined by the local normalization 
condition yielding 

Z = J x{r])e~ { ^ +K+a)v dri. (53) 

We note that Z is a function of 

* = P^-+rt + a, (54) 

expressed by a sort of generalized Laplace transform (this is not exactly a Laplace transform since the variable 77 can 
take positive and negative values): 

r+00 

Z(*)= / x(T?)e-"*tfrj = *(*). (55) 



The coarse-grained angular momentum is given by 

Z(*) 



a = "vvTtTT / VX(v)e n9 dr]. (56) 



It is straightforward to establish that 



<M " Z =F(*) and a 2 = ^=-F'(n (57) 



where 



02 = {a- of = a 2 - g 1 , (58) 

is the centered variance of the local angular momentum distribution (we have noted a n — J prf^drf). According to 
Eq. (fBTl bL F is monotonically decreasing. Therefore, relation l|57la) can be inverted and we get 

piP = -pa--?, (59) 
&7T + pd + a = F- 1 (a) . (60) 

Comparing with Eq. we check explicitly that the coarse-grained flow is a stationary solution of the axisymmetric 
Euler equations. Therefore, for given initial conditions, the statistical theory selects a particular stationary solution 
among all possible ones. We remark that the differential equation for ip, arising from Eq. I|60() and J3J-C involves the 
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inverse F^ 1 of the function determined by Eq. JSTJ-a) while in pure 2D turbulence it involves the direct function F, 
i.e. —Atp = F{j3ip + a). This "inversion" is another striking particularity of our system. 

Comparing Eq. I(6()|l with Eq. 1)24(1 . we note that the equilibrium coarse-grained angular velocity maximizes a 
certain H-function where C is related to F by = —C'(W), i.e. 

C{a) = -J F- 1 (x)dx. (61) 

The H-function S\a] = — f C(a)dr selected by the statistical theory can be viewed as a "generalized entropy" in the 
reduced cf-space |2ft l2ll l22t l23l| . It depends on the initial conditions through the function x(v) which must be related 
to the fine-grained moments of angular momentum (Casimirs). Therefore, in this approach where the constraints 
associated with the Casimirs are treated microcanonically, the generalized entropy in cf-space can only be obtained a 
posteriori, after having solved the full equilibrium equations and related the Lagrange multipliers to the constraints. 
Using (T 2 = — da/d^ according to (f^Tjl -b and = —C'(a) according to f2"4"|) -b. we get the identity 

^ = cWY (62) 

Therefore, at equilibrium, there is a functional relation between the variance o~2 of the distribution and the coarse- 
grained angular momentum cf through the second derivative of the function C. This is similar to the "fluctuation- 
dissipation" theorem |36| . Finally, we note that the most probable value (a)(y,z) of the distribution p{y,z,ri) is such 
that T(vi) — —^77 + lnx(?7) is maximum yielding (ln%)" < and 

( C r) = [(ln X )r 1 W = G(*), (63) 

where G is monotonically decreasing. In general, the most probable value (a) of the distribution l|51|l does not coincide 
with the mean value a. The condition (er) = a is equivalent to 

-(ln^Klnx)']- 1 - (64) 

This equality holds if x is Gaussian. Furthermore, we show in Appendix IH1 that (a) is a stationary solution of the 
axisymmetric Euler equations only when (er) = a. 



C. Particular cases 

Some particular cases of F(^) relationships have been collected in pfT. I2H l2flj. We shall specify different forms of 
function x(v) an d determine the corresponding F(^) and S[W] from the preceding relations. We refer to [2(J for more 
details. In the two-levels case where r] = Xq, Ai, we get 

^-^) = A "+ 1+ A p V 1 A Ao), - ( 65 ) 

In the present case, we need to invert this relation to express £ as a function of cf, hence tp. As discussed above, this 
situation is reversed with respect to pure 2D plane flows. We thus obtain 

^(Ai-Ao) \-pip - 7 - A /x/ 2 fi z y 2p z y \i 



In that case, the generalized entropy in er-space has the form 

Sffl = - [[p\np+(l-p)ln(l-p)}dr, (67) 



with a = p\\ + (1 — p)Xq. This is similar to the Fermi-Dirac entropy. In this two-levels case, the generalized entropy 
S\a] defined by Eqs. (|22(l - ((61|l coincides with the mixing entropy S[p) defined by Eq. (|45|l . This is the only situation 
where we have this equivalence. Taking Xq = and considering the dilute limit a <g; Ai, we get 



Aie- Al * (68) 
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leading to 



-A*V = £ = --f-ln 
Aim 



/iAi / 2/i 2 y 2/i 2 y /i' 



In that case, the generalized entropy in er-space is similar to the Boltzmann entropy 

AT 111 AT' 



SM = - I f Inf dr. 



(69) 



(70) 



II xiv) is a Gaussian, then 

ct = -ct 2 *, (71) 

where the centered variance 02 is a constant. Inverting this relation, we get 

. i . (3 , 7 1 (3 2 lb 8<y a 
[i z a 2 ^<J 2 2 [i z y 2^y /x 

This is the type of mean-field equations that we have considered in Sec. Ill Gl In that case, the generalized entropy in 
cr-space is 



S[a] 



1 



a dr. 



(73) 



which is similar to minus the enstrophy in pure 2D hydrodynamics. 
If x(v) is a decentered Gamma distribution H3-[22fl, then 

cr 2 * 



I + A0-2* ' 



(74) 



1 /2 

where the centered variance 02 is a constant and 2Acr 2 is equal to the skewness of xiv)- Inverting this relation, we 
get 



/3<tp + 7 



1 (3 2 ib f3j a 



(T 2 /x(l - Xpip - A7) 2 [i 2 y 2^ 2 y \i 
In that case, the generalized entropy in Tr-space is 



^1 = -3S 



a— — ln(l + Act) 
A 



dr. 



(75) 



(76) 



D. Relaxation towards the statistical equilibrium state 



We would like now to construct a system of relaxation equations which conserve all the invariants of the inviscid 
dynamics (robust and fragile) and relax towards the statistical equilibrium state. These equations can be used as a 
numerical algorithm to construct the statistical equilibrium state. They also provide a subgrid scale parameterization 
of axisymmetric turbulence. In that context, they can describe the dynamical evolution of the flow on the coarse- 
grained scale. Note that in the coarse-grained formulation, the inviscid approximation is easier to justify, since 
viscosity only acts at very small scales. Following the approach of Robert & Sommeria |24j, these relaxation equations 
can be obtained by using a Maximum Entropy Production Principle (MEPP). 

The equations of evolution for the coarse-grained fields are given by 

-V • 3 a , (77) 

where 3 a and are currents which contain all the information coming from interaction with sub-grid scales. Note 
that we have kept the advective terms because these equations are expected to describe the relaxation of the flow (on 



da 
dt 

dt 



u • Vct = 



u- V£ 
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the coarse-grained scale) towards statistical equilibrium; they are not only numerical algorithms. If we want to keep 
track of the conservation of all the Casimirs (or equivalcntly of the total area of each level of angular momentum) , 
we need to introduce equations of conservation for the density probability p(r, 77, t) of angular momentum. We write 
them as 



-v J, 



(78) 



where J(r, rj,t) is the current of the level 77 of angular momentum. Multiplying Eq. (|78|) by 77 and integrating over 
all the levels, we recover Eq. H77t a) with J CT = J J 77 drj. Furthermore, the conservation of the local normalization 
J pdr\ — 1 imposes 



The time variations of S[p] are given by 



J drj = 0. 



J • V Inp drdr], 



(79) 



(80) 



while the time variations of E and H have been given previously in Eqs. (|27|l - l|28() . Following the MEPP, we maximize 
S with E = H = 0, the normalization constraint 179f) and the additional constraints 



J L dv < C (r,t), -±< C( (r,t). 
Writing the variational principle in the form 



5S - I35E - p5H 
we obtain the optimal currents 



v>'( Y p ) llr,ll i 



x's 



( JJ 



dr = 0, 



J = —D 



Vp + \pp{i] - tf)V I - ) + ppij] - a) V£ 



J 5 = -D'(pVip + iNa), 



(81) 

(82) 

(83) 
(84) 



where we have used Eq. 1)79(1 to obtain the final expression of the current ((83(1 . The current of angular momentum is 
therefore given by 



-D 



(85) 



where 02 is defined in Eq. 1(58(1 . We note that (3{t) and p(t) are time dependent Lagrange multipliers that evolve in 
order to conserve energy and helicity. Their explicit expression is obtained by inserting Eqs. ((84(1 and ((85(1 in the 
constraints l(27|) - (|2^|) . 

We now show that the relaxation equations (JTTJ-b) and 1(78(1 with the currents 1(84(1 and ((83(1 increase the mixing 
entropy ((45(1 until the Gibbs state is reached (ff-theorem). We can write the rate of entropy production l(8U(l in the 
form 



pp{r) - o-)V£ 



drdr] 



(86) 



1 



P / (V - o 7 ) J • V - drdrj + p / (77 - ct) J • V£ drdrj 



Using Eqs. $Q and JHSJl this can be rewritten 



dr + p / J (j ' V£ dr. 



(87) 
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Using Eqs. l(27 )l -l(25 )l . we get 



hence 



S ' ~D~p dV ' J? ' + ^ VZF ) dr ' 



(88) 



(89) 



We conclude that S > provided that 13, 13' > 0. At equilibrium S = yielding J = = 0. Equations l|84|l and 
|J5JJ imply 



Vlnp+i/3(rj-a)V^ ) +/ ,( V -7f)V£ 0. 



(90) 



/3V^ + iiVa = 0. 



(91) 



The second equation is equivalent to 



-4> 



V 

1 1' 



(92) 



On the other hand, for any reference level r/ , Eq. I|9(J|I yields 



Subtracting Eqs. I|90[l - H93|) and integrating, we obtain 



In I — I +^/3(?7-%)-+M(??-%)e = ^), 
VPo/ 2 y 



(93) 



(94) 



which can be written 



P = 



Z(r) 



(95) 



Thus, the stationary solution of the relaxation equations is the Gibbs state (|51|l . 

The relaxation equations are relatively complicated to solve, because we need to solve N coupled PDE, one for each 
level. Alternatively, we can write down a hierarchy of equations for the moments of angular momentum a n . The first 
moment equations of the hierarchy can be written 



Of 



u- V£ 



D 

dz\4y 2 



(96) 
(97) 



We are now led to a complicated closure problem because each equation of the hierarchy involves the next order 
moments. For example, the equation for a involves 02 etc. In the two levels approximation, one has 02 — (cf — 
-V))(-^i — o 7 )- On the other hand, a Gaussian distribution of fluctuations at equilibrium can be obtained by imposing 
that (72 is constant. In these two particular cases, the equations (ffi))l are closed. More generally, we must write 
down the higher moments of the hierarchy and close them with a local maximum entropy principle as proposed by 
Robert & Rosier [37T |. If we implement this procedure up to second moment, it leads to a Gaussian distribution. Its 
implementation to higher moments is difficult. Furthermore, its physical justification is unclear. In practice, we must 
come back to the N coupled PDE for the levels. 
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E. Prior distribution and generalized entropy 



In the statistical approach presented previously, we have assumed that the system is rigorously described by the 
axisymmetric Euler equations so that the conservation of all the Casimirs I n must be taken into account. This 
corresponds to a freely evolving situation. Alternatively, in the case of flows that are forced at small-scales, Ellis 
et al. 28] have proposed to replace the conservation of all the Casimirs by the specification of a prior distribution 
x(jl) encoding the small-scale forcing. This approach has been further developed in Chavanis pl|. In this approach, 
the constraints associated with the (fragile) moments I n >\ are treated canonically instead of microcanonically. By 
contrast, the robust constraints (energy, circulation, helicity, angular momentum) are still treated microcanonically. 
If we view the levels rj of angular momentum as different species of particles, this approach amounts to fixing the 
chemical potentials instead of the total number of particles in each species. The idea is that the ambient medium 
behaves as a reservoir of angular momentum: the small-scale forcing and dissipation affect the conservation of the 
moments of angular momentum I n> \ while fixing instead the canonical variables a n . 

Therefore, in the present situation, the relevant entropy S x is obtained from the mixing entropy (|45(l by making a 
Legendre transform on the fragile moments I n> \ plj . If we assume that the a n in the variational principle H49|) are 
fixed by the "reservoir" (ambient medium), we can define a relative entropy by 



S x = S - 2J at n I n = S 



n>l 



n>l 



a„ / pr/ n dydzdrj. 



(98) 



This is similar to the passage from the entropy S (microcanonical description) to the grand potential f2 = S — p,N 
(grand microcanonical description) in usual thermodynamics when the chemical potential is fixed instead of the 
particle number. Using Eq. (|45|) . we get 



S x = - / p a n Tl 

J L n>l 



dydzdrj. 



Introducing the function (|52ll . we obtain 



(S V 



pin 



x(v) 



dydzdrj. 



(99) 



(100) 



The function x(v) is interpreted as a prior distribution of angular momentum. It is a global distribution of angular 
momentum fixed by the small-scale forcing. It must be regarded as given. In this approach, the statistical equilibrium 
state is obtained by maximizing the relative entropy l|100|) while conserving only the robust constraints. Thus, we 
write the variational problem as 



5S X - f3SE - jiSH - jST - aSI - J ((y, z)s(^J pdr^jdydz = 0. 



(101) 



We can now repeat the calculations of Sec. IIIIBI with almost no modification. The only difference is that we regard 
the a n 's as given. Therefore, the Gibbs state is determined by Eq. (|51f) where x(v) is fixed a priori by the small-scale 
forcing. Recall that in the previous approach (freely evolving flows), it had to be determined a posteriori from the initial 
conditions (assumed known) by a complicated procedure. Here, the specification of xiv) automatically determines 
the function F by Eq. l(57|l and then C by Eq. JJJTJ. Thus, the generalized entropy in cr-space S[a] = — f C(W) dr is 
now determined by the small-scale forcing through the prior xiv) while in the preceding approach it was determined 
by the initial conditions through the Casimirs I n >i- Explicitly, the generalized entropy is expressed as a function of 
X by the formula |23| ]: 



C(a) = -I [(lnxYrH-^dx. (102) 

The equilibrium coarse-grained angular momentum a maximizes the generalized entropy l|22|l - (|102|l at fixed robust 
constraints E, T, H and /. In the present context, the relaxation equations introduced in Sec. Ill Fl can describe 
the relaxation of the coarse-grained flow towards statistical equilibrium once the small-scale forcing has established 
a permanent regime characterized by a prior distribution x(v) determining a generalized entropy S[a]. Since we are 
now interested by the route to equilibrium wc need to restore the advective terms, so we write: 



8a 
dt 



u ■ Vct = v ■ Id 



to+ c4) v (!) 



7f 2 



d_ 

dz \ 4y 2 



V • 



C"(a) 



(103) 
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The physical interpretation of these equations is quite different from (|33|> - In Sec. Ill Fl the relaxation equations 1331) 
provide a numerical algorithm to determine any nonlinearly dynamically stable stationary solution of the axisymmetric 
equations, specified by the convex function C . In this context, only the stationary solution for t — ► +00 matters and 
the evolution towards that state has no physical meaning (it is just the engine of the algorithm). In the present 
section, the relaxation equations (|103fl provide a description of the evolution of the coarse-grained field, for all time 
t, in a medium where a small-scale forcing imposes a prior distribution xiv) ( or a generalized entropy C(<r) in the 
reduced a-space). These equations conserve only the robust constraints and satisfy a generalized maximum entropy 
production principle for the functional S\a]. Finally, the relaxation equations of Sec. IIII Dl provide a description 
of the evolution of the coarse-grained field, for all time t, of a freely evolving system. These equations conserve all 
the constraints (including the Casimirs) and satisfy a maximum entropy production principle for the functional S[p\. 
Note that the relaxation equations I|1U3|) can also be obtained from the moments equations i|9tj|H97|) °f the ordinary 
statistical theory by using the relation I|(j2fl to express 02 as a function of a. In a sense, this relation can be seen as 
a closure relation imposed by a small-scale forcing. Thus, the relaxation equations l|103|) are not simply numerical 
algorithms; they can also provide a parameterization of axisymmetric flows with a small-scale forcing. Their interest 
as numerical algorithms (in the sense of Sec. EFj remains however important in case of incomplete relaxation to 
construct stable stationary solutions of the Euler equation which are not consistent with the statistical theory (in 
cases where the evolution is non-ergodic) in order to reproduce the observations, as discussed in [20l . l2il.l22L I23I ] . 

IV. SUMMARY 

In this paper, we have developed new variational principles to study the structure and the stability of equilibrium 
axisymmetric flows. We have completely characterized the steady states of the inviscid dynamics and found that there 
is an infinite number of solutions. We have shown that each of these steady states extremizes a certain functional 
and that maxima or minima of this functional correspond to nonlinearly dynamically stable states. We have given 
analytical solutions in some simple cases to illustrate our formalism. One of these steady states (non-universal) will 
be reached on the coarse-grained scale as a result of violent relaxation (chaotic mixing) . Our general approach must 
be contrasted from that of other authors who obtained particular solutions of the Navier-Stokes equation by means 
of phenomenological principles (minimum of enstrophy for example in 2D turbulence, Beltramization for MHD and 
axisymmetric flows,...). Such solutions are recovered as particular cases of our formalism but many other solutions 
can emerge in practice depending on the initial conditions, on the route to equilibrium (ergodicity) and on the 
type of forcing. This is why we try to remain very general. Our point is that there is no clear universality in 2D 
or axisymmetric turbulence [2(1 |21L I22L |23| . In a second part, we have developed a thermodynamical approach to 
determine the statistical equilibrium states at some fixed coarse-grained scale. We found that the resulting distribution 
can be divided in two parts: one universal part, coming from the robust constraints, and one non-universal, which 
depends on the initial conditions (Casimirs) for freely evolving systems or on a prior distribution encoding non-ideal 
effects such as forcing and dissipation. Finally, we have derived relaxation equations which can be used either as 
numerical algorithm to compute stable stationary solutions of the axisymmetric Euler equations, or to describe the 
dynamics of the system (freely evolving or forced) at the coarse-grained level. 

The main question regarding the application of our results to realistic systems (such as the ones mentioned in the 
Introduction) is the relevance of the use of the ideal (Euler) equation instead of the true dissipative system. In fact, 
the presence of a small viscosity does not preclude the applicability of our results. First of all, since viscosity acts 
at small scales, its main effect is to erase the fluctuations around the coarse-grained field. Thus, it gives a physical 
support for selecting the coarse-grained field which is at large scales and which is relatively robust against viscosity. 
On the other hand, we have shown that viscosity and coarse-graining act in a similar manner so that they are not in 
opposition. In a very turbulent flow, the diffusion acts only at small scales by dissipating energy. By disregarding the 
details of the fine-grained dynamics, we have a similar process where energy is lost in the small-scales but accumulates 
in the large-scales. On the other hand, we have shown that a (generalized) selective decay principle can be motivated 
either by viscous effects or by coarse-graining. Indeed, a small viscosity or a coarse-graining tend to increase the value 
of the ii-functions (fragile constraints) with only weak modification on the energy, angular momentum, circulation 
and helicity (robust constraints). Therefore, viscous effects do not break the nonlinear dynamical stability results. 
On the contrary, they can precisely explain (together with coarse-graining) how the system can reach a maximum 
of an if -function at fixed robust constraints. Without dissipation (viscosity or coarse-graining) this is not possible 
since the Casimir functionals S[o~] are rigorously conserved by the Euler equation. We believe, however, that the 
main increase of the fragile constraints (like enstrophy) is due to coarse-graining (18| rather than molecular viscosity 
(in classical works on 2D turbulence, it is argued instead that enstrophy is dissipated essentially by viscosity). The 
main difference between viscous and inviscid flows is that inviscid flows tend to a strict stationary solution of the 
Euler equation (on the coarse-grained scale) while, in the presence of a small viscosity, this large-scale structure slowly 
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diffuses and ultimately disapears. However, if v — > 0, this happens on a long time scale that is not of most physical 
interest. Note finally that forcing can act against viscosity and maintain a steady state as for an inviscid evolution. 

The other effect of viscosity, now regarding the statistical mechanics approach, is to break the conservation of 
the Casimirs. This is a problem for the original approach (Sec. IIII B|l where it is assumed that all the Casimirs are 
conserved. However, in the point of view developed in Sec. IIII El we have replaced the specification of the Casimirs by 
a prior distribution of angular momentum. It corresponds to the non-universal part of the distribution of fluctuations 
given by the Gibbs state Q51f) . We have argued that this prior is precisely determined by non- ideal effects such as 
viscosity and forcing (in addition to the initial conditions and the boundary conditions), i.e. by all the complicated 
features of turbulence. Therefore, in this point of view, the existence of a viscosity and a forcing can be taken into 
account phenomenologically in the theory. On the other hand, it should be noted that the effect of coarse-graining is 
similar to a turbulent viscosity. This is best seen in the relaxation equations 111031) which involve a diffusion term with 
a "turbulent viscosity" D. However, our approach shows that the adjunction of a turbulent viscosity to the Eulcr 
equations in order to model turbulence is not sufficient as it breaks the conservation of energy. Therefore, additional 
drift terms arise in the relaxation equations to act against diffusion and lead to a steady state plj . There are other 
pieces of evidence for the claim that the introduction of a coarse- graining procedure is similar to a diffusive process. 
For example, recent numerical simulations have shown that the Euler equation with a high wave-number spectral 
truncation shows similar features as the Navier-Stokes (dissipative) equation |38( . This issue concerning the influence 
of viscosity will be addressed more thoroughly in a second paper, where we confront our prediction to experimental 
data, and use them to derive and characterize the non-universal features of the equilibrium distributions. 

Finally, we will address the changes to be made to account for a global rotation of the system. Taking the rotation 
vector to be aligned in the z-direction, the Coriolis force will only add a term 2Q,u on the left-hand-side of the first 
equation J2J). Then, the conserved quantity will be a' — rv + fir 2 instead of a = rv. This is similar to the use of a 
potential vorticity when doing the statistical mechanics of two-dimensional rotating fluid instead of the usual vorticity 
|36j . Similarly, the right-hand side of the second equation will now be: d z [(v 2 + 2Qvr)/r] — d z [(v 2 + 2Qvr + fl 2 r 2 )/r] = 
d z [a' 2 /r 3 ]. Consequently, all the results in this paper will be valid provided that a' is used instead of a. 
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APPENDIX A: DERIVATION OF CONSERVATION LAWS 

In this Appendix, we prove the conservation laws used in the main text. A cornerstone of the proof is the general 
identity: 

x{<j>, tp} dydz = - J <p{x, ip} dydz , (Al) 

which holds if one of the two fields \ or (p vanishes on the boundary of the domain. 

• Energy conservation: using the equations of motion and assuming that ^ = 0or( = (j = 0on the boundary of 
the domain, we have 

E = J^dydz + J^dydz, (A2) 



1> 



2 V 



£-{i>,<r} dydz 



-£{"0, 0} + ^ {i>, °} - a } j d v dz 
= o, 

where we have used the identity IjAlfl twice to obtain the third line. 

Casimirs conservation: using the equations of motion and ip = or a = on the boundary, we will show that 
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all the moments of a are conserved: 

/da 
a n ~ 1 —dydz, 

= —n J a n ~ l {ip , a} dydz 

= n J ifj{a n ~ 1 ,a} dydz 
= 0, 

where we have used the identity IjAlj! in the third line. 
• Helicity conservation: using the equations of motion, we have 



H 



= ~ / F(a) 



2 V 



J £,F'{a){^,a}dydz 



Then, 



(A3) 



(A4) 



(A5) 



(A6) 



J F(a){^,a}dydz = - J ^{F(a), a} dydz = , 

if a or F(a) vanishes on the boundary of the domain. Therefore, 

H = - J F(a){ip,Z}dydz- J £,F'(a){i/j,a} dydz 

= - J F(a){i;,t}dydz + J £{F(a)^}dydz 

= - J t{F{a)^}dydz + J Z{F(a),iP}dydz 
= 0, 

where we have used identity J8J in the second line and identity (|A1|) in the third line and assumed that F(a) = 
or £ = on the boundary of the domain. 



APPENDIX B: STABILITY OF SOLUTIONS 

In section Til Dl we found that the functions ao(y,z) and £o(y, z ) which extremize the functional H22J) are solutions 
of the following set of equations : 



p-ipo — — £tco — v and 



C'(<to)= 13^ + ^0 + a. 



(Bl) 



However, only maxima of S are nonlinearly dynamically stable. We need therefore to investigate the sign of the second 
order variations of J — S — (3E — \iFl — vY — al . Writing a = o$ + 5a, £ = £o + <5£ and ip = ipo + Sip, one obtains for 
all 5a and 5£: 



5 2 J[a , £o] = -J dydz (c" '(a ) + (5a) 2 + 2»5aSZ + pSifiSt 



(B2) 



with 5£ = —A^Sip. Using the operators curl and Curl defined in [23, one can easily show that the last term can be 
rewritten : 



P / 5ip5£, dydz = p 



Curl 



6if> 



dydz . 



(B3) 
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Putting 5a = in equation {B2j, the condition S 2 J[ao, Co] < thus implies that j3 must be positive. This is at variance 
with pure 2D hydrodynamics, where stable structures can exist at negative temperature and are the most relevant. 
Also, assuming (5£ = 0, one finds that a maximum of J should satisfy the following condition: 



C'M + ^l (<W <u,<h • 1 1 . 



(B4) 



which is trivially fulfilled because C is a convex function. One cannot find a general condition on the value of p in 
order for J 2 J[cr , £o] to be negative but a sufficient condition can be found by using the fact that in (|B2|I . the last 
term in the integral is everywhere positive. Consequently, a sufficient condition for 8 2 J < is: 



dydz 



/ dydz 



> 



C> ) + 




P 



C"(a ) + 



-st; 



>o, 



(B5) 



for all 8a and St;. This condition can obviously be provided if p — 0. A sufficient condition for <to and £oj solution of 
IjBljl . to be maximum of S is thus p = and (3 > 0. However, this is only a very particular case. 



APPENDIX C: STATION ARITY OF (a) 

The most probable value of the distribution l|51|l can be written 

(a) = G(*) = = G 



(CI) 



showing that (cr) = /(i/') is a function of ip alone. We now write the condition under which (a) is a stationary solution 
of the axisymmetric Euler equations. Comparison between Eqs. <|14fl arid i|tj(J|) shows that the following relation must 
hold : 



[i 2 dip 



Using Eq. (|59|l . this can be rewritten: 



P 2 ^P 2 n 7/3 



2^V- 



(C2) 



(C3) 



where c is an integration constant. If we require that a — (a) on the boundary of the domain [ijj = 0), then, using 
Eq. we get / 2 (0) — c = 7 2 /^ 2 . Substituting in Eq. I|C3|I . we find that G o F^ 1 is the identity so that (cr) = a 

(everywhere). This implies that (cr) is a stationary solution of the axisymmetric equations only when it coincides with 
a . 



APPENDIX D: FLUCTUATIONS OF £ 
1. Generalities 



In this Appendix, we try to develop the statistical mechanics approach in the general case, without ignoring the 
fluctuations of £. Since £ is not conserved by the axisymmetric equations (D£/Dt ^ 0), this may invalidate the use of 
a statistical theory to predict its fluctuations, so that our approach is essentially phenomenological and explanatory. 
We introduce p(r,r), v), the density probability of finding the values a = r\ and £ = v in r at equilibrium. Then, the 
coarse-grained fields are a = j pr/drjdv, £ = j pv dr\dv and the local normalization is f pdr/dis = 1. We introduce the 
mixing entropy 



m—f 



plnp dydzdrjdv. 



(Dl) 
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As usual, the fluctuations of ip will be neglected because it is an integrated quantity of the primitive field £. The 
integral constraints can be re-expressed as 

— If 1 f i? 

E = — / pvtp dydzd-qdv + — / p — dydzdrjdv , (D2) 
2 J 4 J y 

iT n = / /wry™ dydzdrjdv , (D3) 
I n = prf 1 dydzdrjdv, , T = pv dydzdrjdv. (D4) 



The most probable distribution at metaequilibrium is therefore obtained by maximizing the entropy at fixed E, H n 
and /„. We introduce Lagrange multipliers and write the variational principle in the form 

SS -I35E~J2 - Yl a » Sl ^ " / C(l/. z ) s ( [ P d vdA dydz = 0. (D5) 

The variations on p yield the Gibbs state 

(D6) 



i p -l3[u^+^]- a G( V )-^F(r,)u 



Z(y, z) 

where jiF(r\) = J2 n Mn 7 ?™ an d oiG{rj) — J2 n a n r f l ■ The "partition function" is determined by the local normalization 
condition yielding 



drjdv. (D7) 



The coarse-grained fields a and ^ are the averaged values of the distribution (|D6(1 . This approach predicts that the 
distribution of fluctuations of pseudo-vorticity is exponential ~ e Hv: z ^i)f so that, depending on the sign of A, it 
diverges either for v — > +oo or v — > — oo. The problem of the smoothness of the vorticity is still an unresolved issue 
for the Navier-Stokes equation related to the existence and uniqueness of its solutions. In two dimensions, it can be 
shown that it is bounded |39| | whereas in three dimensions, very little is known. In our case, which is intermediate 
between these two cases, we can infer that the vorticity will be bounded as we have seen that in the axisymmetric 
case, the vorticity has to vanish in the long time limit (see section Hi E|l . The range of integration for the variable v 
is thus restricted to finite values: \v\ < A. 

We now look for extremum values of the distribution (|D6|I .To study this problem, we write: 



p{y,z,ri,v) = ^ cxp[-J r (7;, u)] , (D8) 

8n 2 

with T = v\j3ip + pF(rj)} + — + aGM . 

Ay 

We start to search for extrema of the distribution in the interior of the domain ] — £ m , £ m [ and call them (a), (£). One 
can check that they obey the equations: 



drj 
dT_ 

dv 



i/fT) + a G'({*))+v{OF'((*))=0, (D9) 



= M + p J F((a))=Q 



These extremum fields (a) and (£) are stationary states of the Eulcr axisymmetric equation (with families indexed 
through the conservation laws) while the averaged states a and £ are not in general. The stability of these extremum 
states can be found in principle by considering the second variations of T . Here, we prefer to use the following trick. 
We introduce the functions : 

Q(c*, £*) = M + nF(a*) , (D10) 

2y 
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Then, for any <t* and , one has 



(3 



Choosing cr, = (er) and 



Z*{y,z. 



■ exp 



= -Q(<j„£*)(z/ - £») - R(v*,Z*)(v <?*? - H(y - &)(ti - v*)F'(a*) (Dll) 

a [G( V ) - G{a*) - {r, - ct*)G'(<7*)] - ^ [F(r?) - F(ct.) - (77 - ct*)F'(ct»)] . 
(£), we have Q = i? = 0, so that the probability function simply becomes: 

(a)) 2 - M („ - (0)(„ - (<7»F'«<t» - a [G(r?) - G((a» - („ - (<7»G"«<7»] (D12) 



, 



- F [F(r,)-F((a))-(,,-(a))F'(H)] 



Due to the e Aly dependence of the density probability, one can check that the extremal states of the Gibbs probability 
distribution are saddle points (stable in one direction, unstable in the other), except when F is constant and f3ip+p,F = 
(leading to £ = 0), in which case they are stable states for positive temperature (3 > 0. The fields (a) and (£) thus 
are not real extrema of the distribution. Moreover, one sees that when F is constant with j3tp + fiF = and G is 
linear, the probability distribution of r\ is a Gaussian in the variable (77 — (c)) 2 and the probability distribution of v 
is uniform. Therefore, the most probable state (a) coincides with the mean state a. However, this is not the generic 
case. 

We now look at possible extrema on the frontier of the domain of integration, for v = ±A. It is obvious that if it 
exists a physical bound on the vorticity, it must depend on the shape of velocity field: A = A(ct, - 0)- Assuming this 
function to be known, we can write the conditions that must satisfy an extremum (co, Co) of p located on the frontier 
of the integration domain: 



d JL = ^ + aG > {ao)+fl( : F>(a )=Q. : 

Of] (7 , So 2y 

Co - ±A(a»- 



(D13) 



We note that these fields are not stationary states of the Euler axisymmetric equation, contrary to (a) and (£) . To 
decide which of the two couples ((a), (£)) or (ao,£o) ^ s ^ ne mos t probable state of the distribution p, one has to 
compare the value of the function T at these two points: 



H°», Co) = Zo\W + M{^o) - <roF'(*o)}} + a{G(a ) - <r G'(a )} 



(D14) 



and Fi(a), (£)) = a[G({a)) (a)G> \{a))] - pF> {{o)){a){t) . (D15) 

From these expressions, it is not possible to decide which one of these two values is the smallest (corresponding to a 
maximum value for p) in the general case. However, in the special case where G oc a and F oc a, we have 

■7>o,£ ) = £o/^, (D16) 



and T{{a), (£)) - - (W, (D17) 

where we have used Eq. (|D9(l -b to obtain the last equality. Since < A, one obtains the most probable state on the 
boundary of the domain of integration. More generally, since corresponds to a saddle point of T 1 the relevant 

solution to consider should be the solution (o"o, Co) where £ reaches its maximum bound. Therefore, this approach 
suggests that the equilibrium states of axisymmetric flows are those that maximize £ (the toroidal component of the 
vorticity). Since the dissipation of kinetic energy is equal to th e sp ace integral of the squared vorticity (see Sec. 
Ill Ell . our conclusion resembles the assumption made by Malkus |40j |. followed by Howard and Busse [42|. who 
calculated bounds on the kinetic energy dissipation for thermal convection problems by maximizing the dissipation 
on a manifold that includes the solutions of the problem. This principle of maximal dissipation has been extended 
to a purely chaotic system by [43J and in this case too, the observed equilibrium solutions are very close to that 
calculated by a maximization of the dissipation. It is however interesting to notice that in these approaches, the 
maximum dissipation is assumed and the shape of the equilibrium field is derived while in our approach we show that 
the equilibrium state (which maximizes the entropy) is the one with a maximal vorticity field. The drawback of our 
approach, however, is that we do not have an explicit form for the equilibrium solution, unless we know how to derive 
the bound on the vorticity. 
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2. Examples 

A few examples can be given to illustrate the points developed above. For simplicity, let us consider first the case 
with F = 1, where rj and v become independent. In such case, the probability distribution function is : 

p= 1 -,( W+fl )-^- aGW) (Dlg) 
Z 

and the partition function factorizes into Z = Z U Z V with : 

Z v = ( + e~"^ + ^dv = sinh[A(/?V + /*)] . (D19) 

Following the discussion of the previous section, we have introduced a symmetrical cut-off A. This situation is similar 
to the Turkington model in 2D turbulence, see [2(J. By integration, one finds 

« - k C ver "^" iv - -rao " + "» ' (D20) 

where 

L(x) = tanh" 1 ^) - - , (D21) 



is the Langevin function. For the other part of Z we get : 

/ + °° „ 2 

e-^r-^Wrfr, . (D22) 
-oo 

• Case G{rf) = r], F = 1: in this case, the extremal state is ip — — /u//3, (a) = —2ay/(3 or cto = —2ay/(3, £o = ±A. 
To derive the mean state, we first compute : 

Zv = 2 (jf) 1 ' (D23) 

The mean state may then be found from : 

<r=— I T)e ^ ar >dri = a- JL = - ~tV ■ (D24) 

Z v J_ x da (j 

The mean state a therefore coincides with the extremal state (a) if tp = — /x//3 or er - The mean state £ is equal to 
zero if ^ = — /i//3 and is lower than |£o| = A (in absolute value) otherwise. 

• Case G(rf) = r\ + fc?? 2 , F = 1: in this case, the extremal state is ip = —/J.//3, (a) = —2ay/((3 + Akay) or 



Co = —2ay/(P + Akay), £o = ±A. For the partition function, we have : 

s (3 + Aaky 



1/2 2 

z v = 2 [ ~r r ) e^yfe , (D25) 



and 

<91nZ„ 2/3y(a - fc) + Aaky 2 {a - 2k) 
da~ ~ ((3 + Aaky) 2 



(D26) 



Therefore, the mean state does not coincide with the extremal state. 

• Case F(rj) — G(rj) = rj: in this case, the extremal state is ip — — (fi/p3)(a), (£) = —ajpk — (/3/2yfj,)(a) or 
£o = — cn/fjt — (/3/2y/i)(To, £o = iA. Integrating the partition function first with respect to rj, we get : 

Z= \j — J dvexp\-vPip-\ - j . (D27) 
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Using this, we find the mean state as : 

r+A 



and 



- 1 Any 



dvv exp 



(a + \xvfy 



1 / 47ry (3 
Z 



+A 

-A 



P 2p </y lj_ x (3 L (3 



]} 



1 I Airy (3 



Z\ (3 fy 



exp 



(« 2 + M 2 A 2 )y 







sinh 



2a[i\y 



+ 



(3 2 j> _ a 
2p?y n ' 



dlnZ 2y. -. 



(D28) 



(D29) 



When the cut-off A is taken into account, the mean state does not coincide with the extremal state. 
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